Two-dimensional Delaunay Triangulations and Voronoi Tessellations in Julia
Department of Mathematics, Imperial College London
2024-10-09
using DelaunayTriangulation, CairoMakie
boundary = [(0.0, 0.0), (0.1, -0.1), (0.2, -0.2), (0.2, 0.0),
(0.8, 0.0), (0.8, -0.2), (0.9, -0.1), (1.0, 0.0),
(0.9, 0.2), (0.7, 0.3), (0.6, 0.2), (0.5, 0.1),
(0.4, 0.2), (0.3, 0.3), (0.1, 0.2), (0.0, 0.0)]
points = [(rand(), 0.5rand() - 0.2) for _ in 1:1000]
boundary_nodes, points = convert_boundary_points_to_indices(
boundary; existing_points=points)
tri = triangulate(points; boundary_nodes)
triplot(tri)
using DelaunayTriangulation, CairoMakie
section_1 = [(0.0, 0.0), (0.1, -0.1)]
section_2 = [(0.1, -0.1), (0.2, -0.2), (0.2, 0.0),
(0.8, 0.0), (0.8, -0.2)]
section_3 = [(0.8, -0.2), (0.9, -0.1), (1.0, 0.0)]
section_4 = [(1.0, 0.0), (0.9, 0.2), (0.7, 0.3), (0.6, 0.2), (0.5, 0.1),
(0.4, 0.2), (0.3, 0.3), (0.1, 0.2), (0.0, 0.0)]
boundary = [section_1, section_2, section_3, section_4]
points = [(rand(), 0.5rand() - 0.2) for _ in 1:1000]
boundary_nodes, points = convert_boundary_points_to_indices(
boundary; existing_points=points)
tri = triangulate(points; boundary_nodes)4-element Vector{Vector{Int64}}:
[1001, 1002]
[1002, 1003, 1004, 1005, 1006]
[1006, 1007, 1008]
[1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1001]
4-element Vector{Vector{Int64}}:
[1001, 1002]
[1002, 1003, 1004, 1005, 1006]
[1006, 1007, 1008]
[1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1001]
Set{Tuple{Int64, Int64}} with 4 elements:
(1004, 1003)
(1006, 1005)
(1003, 1002)
(1005, 1004)
using DelaunayTriangulation
using CairoMakie
θ = LinRange(0, 2pi, 100);
θ = [θ[1:99]; θ[1]];
θr = reverse(θ)
otrx, otry = 10cos.(θ), 10sin.(θ)
innx, inny = 5cos.(θr), 5sin.(θr)
inn2x, inn2y = 2.5cos.(θ), 2.5sin.(θ)
x = [[otrx], [innx], [inn2x]]
y = [[otry], [inny], [inn2y]]
points = [(20rand() - 10, 20rand() - 10)
for _ in 1:1000]
boundary_nodes, points =
convert_boundary_points_to_indices(
x, y; existing_points=points)
tri = triangulate(points;
boundary_nodes)
triplot(tri)
Assess triangles using the radius-edge ratio \(\rho = R/\ell_{\min}\).
\(1/\sqrt 3 \leq \rho <\infty\), with \(1/\sqrt 3\) attained for equilateral triangles.
Related to the smallest angle by
\[ \rho = \frac{1}{2\sin \theta_{\min}}. \]
Thus, controlling \(\rho\) controls the smallest angle.
function refine(tri)
split_all_encroached_segments!(tri)
while has_bad_triangles(tri)
T = get_bad_triangle(tri)
insert_circumcenter!(tri, T)
if circumcenter_encroaches(tri)
undo_insertion!(tri, T)
split_all_encroached_segments!(tri)
else
for V in new_triangles(tri)
ρ, A = radius_edge_ratio(V), area(V)
(ρ > ρₘₐₓ || A > Aₘₐₓ) && mark_bad!(tri, V)
end
end
end
endusing DelaunayTriangulation, CairoMakie
outer = [(0.0, 0.0), (0.5, 0.2499), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]
inner = [(0.25, 0.25), (0.25, 0.75), (0.75, 0.75), (0.75, 0.25), (0.25, 0.25)]
boundary = [[outer], [inner]]
boundary_nodes, points = convert_boundary_points_to_indices(boundary)
segments = Set([(1, 5), (8, 3)])
tri = triangulate(points; boundary_nodes, segments)
refine!(tri; max_area=0.001)
triplot(tri)
Random.seed!(123)
using DelaunayTriangulation, CairoMakie
r₁, r₂, r₃ = 10.0, 5.0, 2.5
c₁ = CircularArc((r₁, 0.0), (r₁, 0.0), (0.0, 0.0))
c₂ = CircularArc((r₂, 0.0), (r₂, 0.0), (0.0, 0.0), positive=false)
c₃ = CircularArc((r₃, 0.0), (r₃, 0.0), (0.0, 0.0))
tri = triangulate(NTuple{2,Float64}[]; boundary_nodes=[[[c₁]], [[c₂]], [[c₃]]])
refine!(tri; max_area=0.001get_area(tri))
triplot(tri)
custom_constraint.domain = [A, B, C, D, A]
river = [A, N, O, P, Q, R, S, T, U, V, W, Z, C, M, L, K, J, I, H, G, F, E, A]
A1, B1 = (50.0, 20.0), (50.0, 25.0)
cir = CircularArc(B1, B1, A1, positive=false)
boundary_nodes, points = convert_boundary_points_to_indices(domain)
river_spl = CatmullRomSpline(river)
tri = triangulate(points; boundary_nodes=[[[cir]], [boundary_nodes]])
Ar = 100 * 30
custom_constraint = (tri, T) -> begin
p, q, r = get_point(tri, triangle_vertices(T)...)
c = (p .+ q .+ r) ./ 3
TA = triangle_area(p, q, r)
is_left(point_position_relative_to_curve(river_spl, c)) && return TA > 1e-4Ar
abs(dist(c, A1) - dist(A1, B1)) / dist(A1, B1) < 1e-1 && return TA > 1e-4Ar
return false
end
refine!(tri; custom_constraint=custom_constraint, max_area=1e-1Ar)
triplot(tri)